function [U,W,J,meanprices,Gdest,Jf,Gunc] = alg_1(max_iter,cs,cu,c_inv,sigma,M,beta,ksi,g_euler,lambda,lambda_f,gamma_e,r,delta,Uold,K,w_ii,sigma_exp,hq)

difU = 100;
epsU = 10^-6;
iter = 0;
W = zeros(M,M);

while difU > epsU && iter<max_iter    
     %-- Value of a traveling ship
   
    for i = 1 : M        
        for j = 1 : M            
            if i~=j
                W(i,j) = (1/(1-beta*(1-ksi(i,j))))*(-cs(i,j) + ksi(i,j)*beta*Uold(j));
            else
                W(i,j) = beta*Uold(i);
            end
        end
    end    
    %-- Value of an unmatched ship 
    J = (log(  sum( exp(W ),2 )) + g_euler);
    
    if sum(J)==-Inf || sum(J)==Inf
        a = mean(W,2);  % if it doesn't work, try taking out the min
        a = repmat(a,[1,M]);
        J = (a(:,1) + log(  sum( exp(W-a),2 )) + g_euler);       
    end
    
    if sum(J)==-Inf
        a = min(W,[],2);  % if it doesn't work, try taking out the min
        a = repmat(a,[1,M]);
        J = (a(:,1) + log(  sum( exp(W-a),2 )) + g_euler);       
    end
    
    if sum(J)==Inf
       a = max(W,[],2);  % if it doesn't work, try taking out the min
        a = repmat(a,[1,M]);
        J = (a(:,1) + log(  sum( exp(W-a),2 )) + g_euler); 
    end
    
    rq = r-hq;
    %-- Update the expected prices by origin and destination (see Eq. 10 in the paper)
    
    meanprices = zeros(M,M);
    
    for i = 1:M
        t1 = (1-gamma_e(i))./(1-beta*delta*(1-gamma_e(i)*lambda_f(i)));
        t1 = t1.*(beta.*delta.*c_inv(i,:));
        
        t2 = gamma_e(i).*(1-beta.*delta.*(1-lambda_f(i)))./(1-beta*delta*(1-gamma_e(i)*lambda_f(i)));
        t2 = t2.*(J(i)-W(i,:))./sigma;
        
        t3 = (1-gamma_e(i))./(1-beta*delta*(1-gamma_e(i)*lambda_f(i)));
        t3 = t3.*(1-beta*delta).*rq(i,:);
        meanprices(i,:) = t1+t2+t3;
    end
    
    meanprices(eye(size(meanprices))==1) = nan;
    
    %-- Compute the exporter's value function, given the new prices
    
    Jf = zeros(M,M);
    
    for i=1:M
        
        alpha = (1 - delta.*beta.*(1-lambda_f(i)));
        
        Jf(i,:) = lambda_f(i).*(rq(i,:)-meanprices(i,:))-c_inv(i,:);
        
        Jf(i,:) = Jf(i,:)./alpha;
        
    end
    
    Jf(eye(size(Jf))==1)=nan;
    %-- Exporters' destination shares
    
    vvs = cat(2,w_ii,Jf-K)./sigma_exp;
    g = exp(vvs);
    Gunc = g./nansum(g,2);
    Gdest = Gunc(:,2:end)./(1-Gunc(:,1));
    %-- Expected prices for the ship, given the new destination shares
    
    Ev_price = nansum(Gdest.*meanprices,2);    
    Ev_price = sigma.*Ev_price;
    
    %-- Value of a ship at the beginning of the period 
    
    Unew = -cu + (1-lambda).*(J) + lambda.*nansum(Gdest.*W,2) + lambda.*Ev_price;  
    
    %-- Check convergence
    
    iter = iter+1;
    difU = max(abs(Uold-Unew));
    
    Uold = Unew;
    
end

U = Unew;

end

